Pricing Step Options under the CEV and other Solvable 

Diffusion Models 

G. Campolieti, R. Makarov, and K. Wouterloot 
Department of Mathematics, Wilfrid Laurier University 
Waterloo, Ontario, Canada 
gcampoliSwlu. ca, rmakarovOwlu . ca and kwouterloot@gmail . com 

January 25, 2013 

Abstract 

We consider a special family of occupation-time derivatives, namely proportional step op- 
tions introduced by Linetsky in [Math. Finance, 9, 55-96 (1999)]. We develop new closed-form 
spectral expansions for pricing such options under a class of nonlinear volatility diffusion pro- 
cesses which includes the constant-elasticity-of-variance (CEV) model as an example. In par- 
ticular, we derive a general analytically exact expression for the resolvent kernel (i.e. Green's 
function) of such processes with killing at an exponential stopping time (independent of the 
process) of occupation above or below a fixed level. Moreover, we succeed in Laplace inverting 
the resolvent kernel and thereby derive newly closed-form spectral expansion formulae for the 
transition probability density of such processes with killing. The spectral expansion formulae 
are rapidly convergent and easy-to-implement as they are based simply on knowledge of a pair 
of fundamental solutions for an underlying solvable diffusion process. We apply the spectral ex- 
pansion formulae to the pricing of proportional step options for four specific families of solvable 
nonlinear diffusion asset price models that include the CEV diffusion model and three other 
multi-parameter state-dependent local volatility confluent hypergeometric diffusion processes. 

1 Introduction 

Consider a continuous-time stochastic asset (e.g. stock) price process S = {St}t>o- We recall that 
the occupation times A^'^ = A^'g of the process S below a given level L > 0, during the time 
interval [0,T], is defined by 

Jo 

and, in a similar way, the occupation time for staying above level L is 

fT 

Jo 

Occupation time derivatives were introduced as a more flexible alternative to standard barrier 
and lookback options. Many types of occupation time options have been proposed such as the step, 
a-quantile, Parisian, and corridor options (e.g. see [SJinillS]). A family of proportional step options 
was introduced in |13) as a flexible alternative to knock-out barrier options. Step options do not 
lose their value when the barrier is reached. Instead, the payoff function of such options depends 
continuously on the time that the underlying price spends below or above a given (barrier) level. 
For a proportional step option, its payoff at maturity is defined as the payoff of a vanilla option 
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discounted by a factor exp (^aAx), where At is an occupation time. For maturity time T > and 
given level L > 0, the payoff functions of the proportional "down-and-out" (/stop) ^-^d "up-and-out" 



Here, a > is a parameter and /{St) is a payoff of a vanilla option, e.g. {St ~ K)^ for a European 
call or {K — St)~^ for a put for a given strike K > 0. We note that for a = the proportional step 
option is simply a vanilla European option with payoff /{St)- For any occupation time At > 0, the 
factor e~"'^^ is non-increasing in a. Recall the maximum AIt = supjS't : < t < T} and minimum 
rriT = inf{5t : < < < T} of the process up to time T. By continuity of the process S on ]R,+ , for 
So < L we have e^"'^^^ > ImtKL and e^"^r'^ _> a.s., as a — ?> oo. Similarly, for > L, 

y>i a.s., as a — og. We hence observe that proportional step 
options may be regarded as an interesting alternative to knock-out barrier options with less risky 
payoffs but approach their knock-out barrier counterparts as the discount penalty factor a — > oo. 

Closed form pricing and hedging formulae have been found for many types of occupation time 
derivatives under the assumption that the underlying asset price follows a geometric Brownian 
motion (e.g. see [3 El [TOl [13] ) . However, analytical pricing of occupation time options is a non- 
trivial problem for nonlinear diffusion asset price models. Recently, a Laplace transform-based 
approach to price occupation time derivatives under Kou's double exponential jump diffusion model 
was presented in [3 . In |12 , the double Laplace transform of the joint probability density function 
(PDF) of the asset value and occupation time was obtained for the CEV diffusion model. 

Li this paper, we present an analytical method for pricing proportional step options under a class 
of solvable diffusion models. The models considered include the constant elasticity of variance (CEV) 
diffusion model and other confluent hypergeometric processes. Our approach uses a closed-form 
spectral expansion for the transition density of the process with killing at an exponential stopping 
time of occupation for the process above or below a fixed level. This is essentially the Feynman-Kac 
theorem combined with an analytical inversion method for the Laplace transform of the Green's 
function. The Feynman-Kac theorem has been used in the past to find an expression for joint 
probability distribution of St and A^'^, assuming that {S'f}t>o follows a Brownian motion. It 
has also been used to generate analytical prices for occupation time options assuming that the 
underlying asset follows a GBM in [TUIIIS], as well as Kou's model in [3]. In [T^], the Feynman- 
Kac approach is proposed for deriving the double Laplace transform of the joint probability density 
function (PDF) of 5*^ and A^'^ , assuming that the underlying asset follows a CEV process. In 
principle, this method can be used to price proportional step options and any other options whose 
payoff functions depend only on St and Aif'^ , and it can easily be generalized to a class of solvable 
diffusions. In this paper, we follow the approach similar to that of [T^. In contrast to [12 , we are 
able to derive computationally tractable pricing formulae for proportional step options. By applying 
the Feynman-Kac formula, we obtain the Laplace transform of the transition PDF of the asset price 
process with killing at an exponential stopping time of occupation for the process above or below 
a fixed level. The residue theorem allows us to invert the Laplace transform. The resulting pricing 
formula is given in the form of an integral of a spectral series expansion. 

Our pricing method for proportional step options does not rely on computing the joint PDF of 
St and A^'^ . Although, such a joint PDF is also obtainable by a single Laplace inversion of the 
above mentioned transition PDF. Moreover, the approach in this paper can be readily generalized 
to other step options with a different structure of the occupation time (e.g. corridor options) as the 
the relevant transition densities follow by the same spectral expansion methodology presented in 
this paper. For example, it can be applied to the case of a double-barrier proportional step option 
where the payoff depends on the occupation time of the underlying process in between two barriers 



(/stop) step options for the asset price process {St}t>o take the form: 




(1) 
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Li and L2 with Li < L2 (see Such an occupation time is given by 

^T^S ^ = / lLi<St<L2dt. 

Jo 

We also note that the spectral expansion approach in this paper is generally applicable to pricing 
step options in the presence of knock-out barriers for the solvable models considered. 

The organization of our paper is as follows. In Section [2] we present the general framework for 
the no-arbitrage evaluation of step options and their deltas under solvable diffusion models. This 
section also contains our main result on the Green's functions and the corresponding closed-form 
spectral expansions for the transition probability densities with killing at an exponential stopping 
time of occupation above or below a fixed level. In Section [3] we give explicit analytical expressions 
for the CEV model [7 and the other solvable nonlinear local volatility models [5]. Section |4] presents 
the methodology for computing the step options. In particular, the actual implementation of the 
spectral expansions is given in Subsection |4.1| and a Monte Garlo bridge sampling approximation 
approach is presented in Subsection |4.2| Numerical results for pricing and hedging proportional 
step options under the CEV model and other classes of confluent hypergeometric diffusion models 
are presented in Section |4.3| and some conclusions are drawn in Section [5j Supporting proofs of our 
main theoretical results are contained in the Appendix. 

2 Analytical Pricing of Step Options 

Consider a filtered probability space (17, J", { J"i}t>o, P) and an asset price process S = {S't}f>o, as 
a non-negative diffusion process started at 6*0 = 5 G 11+ adapted to its natural filtration {J-t\t>o 
and having the infinitesimal generator 

{Qsf){S) ■■= \ <j'{S)f"{S) + (r - q)Sf{S) (2) 

acting on a bounded twice continuously differentiable function / : R-|_ — > 11. Here, r > and q > 
are respectively the constant risk-free interest rate and the constant dividend yield rate. We will 
denote P as an equivalent martingale measure (i.e. risk-neutral probability measure) with discounted 
process {e~^''~'i^*S't}t>o as a P-martingale. In all pricing models considered in this paper this property 
holds true. The asset price obeys the stochastic differential equation dSt = (r — q)dt + (T{St)dWt, 
where {Wt}t>o is a standard P-Brownian motion. In the sequel, we shall simply set the dividend 
yield q = 0. 

The boundary behavior of the S-process at the origin and infinity depends on the growth behavior 
of <j{S) as 5 — >■ and 5* — >■ 00, respectively. For the models considered in this paper, 5* = 00 is 
a natural boundary. The point S — can be either natural, or a regular boundary (which we shall 
specify as killing) or an exit boundary. If the equity price process hits the zero boundary before 
the maturity date, the equity goes to bankruptcy and a derivative on the equity becomes worthless. 
Therefore, for a model with default, the payoff in equation ([I]) takes the form 

/slp = e-''^-:-V(^T)lT<.o, 

where tq is the first hitting time at zero. We note that this is equivalent to defining the payoff by 
equation ([ij where f{d) = and the S-process is given by the cemetery state d upon hitting the 
origin, i.e. St = d for T > tq. 

It is a typical situation when a solvable model can be obtained from another solvable underlying 
process, say a diffusion X, by simply applying a change of variables. Such a transformation simplifies 
our formulation and it proves convenient to work directly with the underlying diffusion X = {Xt}t>o 
instead of the asset price process S. We assume that the asset price process S is given by a strictly 
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increasing smooth C^(I) map F : I — > ]R+ of the process X, i.e. St — f{Xt),t > 0. The process 
X is taken to be a one-dimensional time-homogeneous regular diffusion on an interval I = {l,r), 
— oo < I < r < oo. We denote the inverse map by X = F~^. Since F is strictly increasing, the origin 
S = and the point at infinity S = oo have the same boundary classification as the respective left, 
I, and right, r, endpoints of X. Moreover, the occupation times for processes S and X are simply 
related: 



A 



L.+ 



1 



du ; A 



t,s 



A 



t.x 



1 



where F{£) = L, £ = X(L) for any level £ £ X, i.e. I < £ < r. The payoff in equation ([T]) is 
then equivalently given in terms of the occupation times for X and its terminal value: f^ep — 

exp(— x)^(^t) where h{x) := f{f{x)). We note also that the event corresponding to the S- 
process hitting zero is equivalent to the X-process hitting the left boundary I, at which time the 
process X is sent to the cemetery state. Hence, for an exit or regular killing boundary I we have 
Xt = d for T > tq and we set h = 0. 

Throughout we respectively denote the (risk-neutral) expectation operator and probability mea- 
sure for X started at Xq = a: by E^, and P^;. By the map we have the spot S = Sq = f{x) and 
X = X(5'). The no-arbitrage prices V^^^ of proportional step options (with constant interest rate 
and zero dividend) take the form 



,-rT 



^hiXr) 
h{v)P»^{T;x,y) dy. 



(3) 



Hence, the problem of pricing proportional step options is reduced to evaluating an integral involving 
the transition PDF p^'"*" for the process X additionally killed according to its (a-proportional) 
occupation time above level £ or for killing according to its occupation time below level £. 

Let T ~ Exp(A) be an exponentially distributed stopping time with parameter A > 0, independent 
of the process X. For A > 0, the Green's functions give us expressions for the expectations of 
important related functionals involving the process at the exponentially stopped time r: 



E, 
E, 



■> ; Xr e dy 



= E, 
= E, 



exp 
exp I —a 



a I tx^>idu ) ; Xr edy 



1 



du I ; Xr G dy 



^\G', + {x,y,X)dy, (4) 
= XG'--{x,y,X)dy. (5) 



The transition PDFs, Pa^{t] x, y), x, y G I, i > 0, are given by the Laplace inverse (with respect to 
complex A) of the respective Green's functions: 



Pa^{i'^x,y) dy = E^^ 



.'x -XtCzdy 



= L 



E, 



-X -.Xr&dy {t) 



= Ll\G'^^{x,y,\)){t)dy. 



(6) 



The Green's functions G'„'*(a;, y, A) = e 



{Pa^ {t; X, y)) (A) are given by the 



Laplace transform of the respective transition PDFs. Note that, throughout our paper, the Green's 
functions and transition PDFs are defined with respect to the Lebesgue measure. In what follows we 
firstly derive closed-form analytical expressions for G^^^{x,y,X), and hence for the expectations in 
Q and ([5|. This is the result in Lemma [l] Secondly, we proceed to analytically invert the Laplace 
transform in ^ and obtain closed-form spectral expansions for p^^^{t;x,y) in Proposition [l] We 
note that our results are valid for quite general diffusions. The transition PDFs are later used to 
compute the option prices given by the integral in ([S]). 
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At this point we remark on the important connection between the above transition PDF and the 
joint PDF of the occupation time Af^ G [0,t] and the terminal value Xt e I for the process started 
at Xq — X. In particular, using the joint PDF (for either occupation time above or below level C) 

(AjIx e dw , e Ay^ y\x)duAy 

the expectation in (|6]) takes the equivalent form 

[e-"^*:x ; X, e dy] - [j\-'^^v'txS^A^)^^ dy = £„ (p^J^^,^ (u, j/jx)) (a) dy. 

Hence, p^^j^^ (u, and pl^^{t;x,y) are Laplace transforms of one another. In particular. 

Alternatively, this joint PDF can also be expressed as a double inverse Laplace transform of the 
Green's function: 

/At^^iu,y\x)=C-' [c^'{G'„^{x,y,X)){t)) (u). (8) 

Since we are able to analytically invert the Laplace transform with respect to A, the joint density is 
given by a single Laplace inverse of the transition PDF p^^{t;x,y) as in ([t]). This single Laplace 
inversion operation can be readily performed numerically. Being given the joint PDF pAt,Xt: the 
no-arbitrage price of a general derivative contract whose payoff is a function of the occupation time 
and the terminal asset price, /{At^S: St), can be computed as follows: 

V{S, T) - c-'-^E [/(At,s, 5t) I 5*0 - 5] - c->-^E, [/(At,x, F(Xt))] 

f{u,y)pAT,XT{u,y\x)dudy (9) 



0-'-^ / fiu,y)C-'{p'^^iT;x,y))iu)dudy. 



Here we also note that this double integral reduces to the single integral in ([s]) in the particular 
case of pricing step options. Clearly, equation ([s]) gives the more direct method for computing step 
options, both analytically and numerically. If we only had the double Laplace transform of the PDF 
PAt,Xt, the calculation of the derivative price would involve possibly a two-dimensional integral in 
the double Laplace inverse to obtain pAt,Xt and then the two-dimensional integral over time and 
space. 

We now proceed with the details of the diffusion processes that we consider in this paper. We 
begin by assuming that the underlying diffusion X has the infinitesimal generator 

{gf)ix):=^b'ix)rix)+a{x)nx) (10) 

acting on a bounded twice continuously differentiable function / : I — > R. The (infinitesimal) drift 
and diffusion coefhcient functions are assumed smooth with continuous a{x), a'{x), b{x) > 0, b"(x) on 
the open interval I. The diffusion X hence has smooth scale and speed density functions respectively 
defined by 

,{x) := exp 1 ^dz) and m{x) := (11) 



We recall that the Green's function G{x, y, A) for process X (with generator in ( 10 )) has the standard 
form (e.g. see [2]): 

G(x, y. A) = {Wxr^My)Mx A y)4>x{x V y), (12) 
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X A y := inm{x,y} and xW y :— max{a;,?/}. The pair of functions {ipx,4>x\ solve Qif{x) — \ip(x). 
Their Wronskian is given by W[(j)\,'ili\\{x) — 'W\s[x) with Wx depending only on A. [Throughout, 
the Wronskian of two functions is denoted by := f{x)g'{x) — g{x)f'{x).] These functions 

are uniquely characterized (within a multiplicative constant) by requiring that, for real values of 
A > 0, V'A and 0a are respectively increasing and decreasing functions on I and by additionally posing 
boundary conditions at regular (non-singular) boundaries of X (see |5j). For a regular left boundary 
'>pail+) — iil is specified as killing or -^(j^ '^^"Jx^^ = if Hs specified as reflecting and included 
in the state space. If Hs a singular boundary, the functions have the following boundary properties: if 
I is entrance(=entrance-not-exit), then 'ipa{l+) > 0, s(i+) '^^°dx^'^ = 0; if / is exit(=exit-not-entrance). 



then -00 (^+) = Oi 5(1+) dx > 0; if Z is a natural boundary, then ^pa{l+) = 0, 
Analogous conditions hold for the right boundary: ipa{r~) > 0, 
Mr-) = 0, 



d<tic.(r~) 



,(r-) 



dx 



0. 



dx 



= if r is entrance; 

— dx — < if r is exit; Va(f— ) — 0, s(r-) ^'^"dx"^ = if r is a natural boundary. 
In what follows we shall deal with Wronskians of fundamental solutions with differing values 
of the eigenvalue parameter. Hence, it is convenient to adopt a slightly more compact notation as 
follows. Let A, 7 G C, then we define: 



W^:^{x) ■.^W[cbx,^^]{x), 



W[<l>x,i^j]{x), Wt:^{x) := W[iJx,Mi^), 



(13) 



where W^'f{x) 



7,A V-^; - ~W^;t{x) and for 7 



A we recover the above Wronskian W^''f{x) 



Wxsix). 



Lemma 1. Let X have the Green's function in j^lS^ with tpx and (j)x specified as above. Then, the 
Green's function G ~ G^'+(x, j/, A), solving {Q^ — X)G = —S{x — y) with infinitesimal generator 
:= {Q — alx>e), a G H, and having boundary conditions as in G, is given by: 



G'„+{x,y,\)=l 



G{x,y,\) 



\+a,\ 



Wx 



m{y)'4jxix)(t)x+a{y), 



m{y)(l)x+aix)'ipx{y), 



<f+J^) , ,(/.A+o(a^)</>A+a(t/) 



G{x,y,\ + a) + 



M^aX.aW 



A+a 



x<l,y<l, 



x<l,y>l, 



x>i,y< 



, x>e,y>e. 



(14) 



Silimarly, G — G^' (x, y, A) solving [Q,^ — A)G = —5{x — y) with infinitesimal generator 
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(Q — alx<e), a € M,, and having boundary conditions as in G, is given by: 



s{£) 

... "^(y)V'A+a {x)(l)x {y) , 



Gi' (x, ?/, A) = < 



G{x,y,X) - 



rn(y)</'A(a;)V'A+a(y), 



tn W . 



x<e,y>e, 



x>e,y<e, 



x>£,y> 



(15) 



Proof. See Appendix ] A. 1| 



□ 



Note that in the trivial case where a = we recover the Green's function for process X on I, i.e. 
GQ^{x,y,X) = G{x,y,\). In the hmit a — > oo, it can be shown that the above Green's functions 
respectively recover those for the process killed at a lower or upper level i (see the remark just after 
the proof of Lemma [T]) . 

An analytical inversion of the respective Laplace transforms leads to closed-form spectral ex- 



pansions for the transition PDFs. The Green's functions in equations (14) and (15 1 are functions 



of complex A G C and the respective Laplace inverses are given by the Bromwich contour integral: 
p^^{t;x,y) = J^_^^ e^*G^^^{x,y,\)dX, where all singularities of the respective Green's func- 
tions lie to the left of the Bromwich line RcA = c. The spectral expansion for p^^^{t;x,y) can be 
obtained by appropriately closing the Bromwich contour and applying Cauchy's Residue Theorem. 
The residue contributions from simple poles give rise to the discrete part of the spectral expansion, 
while continuous parts of the spectral expansion can arise as integrals over branch cut discontinuities 
of G in the complex A plane to the left of the Bromwich line. The analytical form of the spectral 
expansion of p^'^ {t; x, y) clearly depends upon the singularity structure of G^'* (x, y, A) . 

If the Green's function is a meromorphic function that is analytic in A with the exception of 
a countable number of isolated simple poles then the corresponding transition PDF has a discrete 
spectral expansion. In this case we are in so-called Spectral Category I. We refer to [T5] for a general 
discussion and summary of the possible spectral categories and their relation to the boundary clas- 
sification of the cndpoints for a one-dimensional time-homogeneous diffusion. The spectral category 
and the analytic properties of the fundamental pair {ipx, 4'\} are intimately related to the boundary 
classification. We remark that the boundary classification (i.e. natural, exit, entrance or regular) 
of the endpoints I, r is the same for process X, with generator Q, as for the processes X^'^ defined 
by the respective generators in Lemma [l] This fact is easily shown since all processes have the 
same above scale and speed measure and the additional killing rate, alx>i or alj-^i, is a piece- 
wise constant function of x, and hence does not affect the Feller conditions. Moreover, the two 
Sturm-Liouville (SL) operators —G^ and —Q differ only by a piecewise constant function atx>e or 
al.j;<e- The so-called potential function in the Liouville normal form of the SL equation associated 
to the respective sets of operators —G^ and also differs only by a piecewise constant function. 
It follows that an endpoint e G {l,r} is non-oscillatory (NONOSC) for the process X if and only if 
it is NONOSC for processes X^'='=. Moreover, an endpoint e G {l,r} that is 0-NO (oscillatory/non- 
oscillatory) for X will also be O-NO for processes X^'* with spectral cutoff that may be shifted by 
the amount a. 

If both endpoints {l,r} are NONOSC, the eigenspectrum of the SL operator is simple (non- 
negative for a > 0) and purely discrete, i.e. we are in Spectral Category I. Non-natural (regular. 
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exit or entrance) boundaries are always NONOSC, while natural boundaries may be NONOSC or 
0-NO with some spectral cutoff value. In what follows we shall assume that Spectral Category I 
holds where the endpoints {I, r} of the diffusion X are NONOSC. The Green's functions G^'^(a;, y, A) 
in Lemma [l] are meromorphic functions of A with simple poles at A = — A„,n = 1,2,..., where the 
set of eigenvalues {A„}„>i corresponding to the SL operator —5+ (or —G^) form a monotonically 
increasing sequence of real values (non- negative for a > 0): Xn oo as n oo. The following 
result gives us the spectral expansions for p^^^{t;x,y) in the general case of Spectral Category I. 

Proposition 1. Assume Spectral Category I with endpoints of the diffusion X as NONOSC. Then, 
the Green functions G^'* (x, A) in (14-) and (15) are meromorphic functions of X and the corre- 
sponding transition PDFs (i.e. Laplace transform inverses) have the discrete spectral expansions 



X, y) = m(y) £ e-^"*0l-± (y). 



(16) 



For p\ 



4>'n%{x)4>li%{y) 



^^V'-A„(a;)(/)_^^_^„(?/), 



7rp'/'_A„+a(2;)V'_Aj2/), 



X < £, y < 

x < y > ^, 
x>l,y<t, 



(17) 



where Ci 



solving W'^f {£) = (). 

Forp^^^: 



^A+Q a(^)Ia=-a with eigenvalues {A„}„>i as the set of increasing simple zeros 



— Xn-\-Oi, — Xn 



'n.Q 



w 



V'-A„+aWV-'-A„+a(2/)' X < £, y < i, 

'?^-A„(^)^-A„+a(y)' x>e,y<e, 



(18) 



(t)_l(x)(j)_i{y), 



x>i,y>i. 



A„ 



where ^ := ^^^a 'a^q(^)Ia=-a "with eigenvalues {A„}„>i as the set of increasing simple 



solving Wf^f ^^(£) = 0. 
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Proof. See Appendix ] A. 2 1 □ 

By using the spectral expansions in Proposition [T] within the integral in ([s]), we are now able to 
compute the no-arbitrage prices of various proportional step options. In this paper, we are interested 
in pricing proportional step options with call and put payoffs, where /{St) is either the vanilla call 
payoff {St — K)^ or the vanilla put payoff {K — St)~^- However, our method can be used for 
any well-behaved payoff function /. From equation (|3|, the respective pricing of the vanilla step 
(up/down) call and put options with level L are given by the integrals 



CLp{S,T) = e-^-^ f (F(y) - K)f^^{T-x,y)Ay, 

J X 



step 



(5,T) 



-rT 



{K-V{y))ff^{T-x,y) Ay 



(19) 
(20) 



where = X{K),x = X{S),l = X(L),; = X(0),r = X(cx)), X = p-i. Note that the deltas (and 
gammas) of step options are also readily computed. For instance, by differentiating equation ([s]) 
with respect to S we obtain the delta: 



step 



d_ 
dS 



F±ep(5, T) = e--^X'(^) /(F(y)) (T; x, y) dy. 



(21) 



As an alternative to spectral expansions, one can make use of the numerical Laplace inversion 
when computing the no-arbitrage prices of step options. Such an approach leads to computing the 
following two-dimensional integral: 



-rT 



-tT 



%)£-i(Gl^±(x,y,A))(T)dy 



h{y) 



— rT ^ cT 

e -e 



1 

Hy) 



c+'ioo 



e^* Gi'±(a;,y,A)dA 



dy (or) 



i{ut)^^^G^^^{x,y,c + iu)^ du 



dy, 



where all singularities are to left of the Bromwich line SRA = c. The numerical inversion of the 
Laplace transform can be done by employing the Euler numerical algorithm. However, the Green's 
functions of models considered in this paper are formed of hypergeometric functions which, for com- 
plex arguments, are quite complicated to compute. This circumstance makes the numerical Laplace 
inversion possible, yet not as simple a computational procedure. The use of spectral expansions, 
which are uniform rapidly convergent series, gives us the added advantage of computing p^'^ to very 
high accuracy and efficiency, especially for larger values of T. 



3 Pricing under Solvable Models 

According to Proposition [l] to apply the spectral expansion method and find no-arbitrage prices 
of step options under a diffusion asset price model (specified by the infinitesimal generator Q and 
defined on state space I), we only need to know two linearly independent fundamental solutions, 
4>\{x) and %l)x{x)^ that solve the equation {Q f){x) — Xf{x), A G C, a; S X, subject to respective 
boundary conditions at the endpoints of I. Other ingredients such as derivatives and Wronskians 
of the fundamental solutions can be computed analytically or numerically via finite differences. 
Since the zeros {A„}„>i w.r.t. A of either Wronskian Wt'x+a -xi^) ^-x-x+ai^)^ which can 
be computed numerically, converge to the zeros {A„}„>i of W-x — W^[0-a, V'-a](2;)/s(x) in the 
limit a — > 0+, we are also interested in the distribution of {A„}„>i. First, these zeros can be used 
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as initial guesses for computing {A,i}„>i. Second, we analyze whether the zeros grow linearly or 
quadratically. In the latter case, the spectral expansion series converges more rapidly than in the 
former case and fewer terms are required to achieve a high accuracy of computations. 

In the sequel, we present four analytically solvable asset price models with state-dependent 
volatility functions. The first model to be considered is the well-known constant elasticity of vari- 
ance (CEV) diffusion model. Three other alternative models are constructed by using the "diffusion 
canonical transformation" (see ^ |HJ [S] for details). In particular, we consider three examples of 
hypergeometric diffusion models namely the Bessel-K, confluent-U, and UOU models respectively 
constructed from the squared Bessel, CIR, and Ornstein-Uhlenbeck processes by using the afore- 
mentioned method. For all four models considered in the sequel, the discounted asset price process 
{e^''*iS't}t>o is a martingale under the risk-neutral probability measure P where we assume zero 
dividend on the stock. 



3.1 The CEV Diffusion Model 

The constant elasticity of variance (CEV) diffusion, S = {S'tjj^g G R+ is defined by the infinitesimal 
generator {g^ f){S) = \ 5'^S^^+^f"{S) + rSf'{S) with 6 > Oand r g R. Here we take r > and 
assume (3 < 0. Hence, the point 5" = oo is a natural boundary. For /? < —1/2, the point 5' = is 
a regular boundary, which we specify as killing, and for — 1/2 < /? < it is an exit boundary. 

Recall that the Cox-IngersoU-Ross (CIR) model (known also as the squared radial Ornstein- 
Uhlenbeck process, see |2]) has the infinitesimal generator 

{G f)i^) = \^^^f"{^) + (70 - lix)f'{x), (22) 

with constant parameters 70, 71, and v > Q. The strictly increasing mapping X(5) = S^^ 13^^ S^^^ 
reduces the CEV process to the CIR model with the parameters z/ = 2, 70 = 2 + ^, and 71 = 2r/3. 

It is convenient to define the parameter /i = ^ ^ ^ ~ ~ ^ ~ respective speed and scale 

densities are m(a;) = ^x'^e~'^^''^/^ and s{x) — x^^^^o^^^l"^ . The left endpoint ^ = is regular killing 
for /i G (—1, 0) and is exit for /i < —1; the right endpoint r = 00 is NONOSC natural. Hence, both 
endpoints are NONOSC and Proposition [l] is applicable. 

The respective fundamental solutions for the CIR process are 



V l7i| 2 
^A(2:)=xl''le^i-/2_^ (^1 + _^,1 + |^|,M^ 



(23) 



where M.{a,h,z) and V({a,b,z) are the confluent hypergeometric functions (see [T]). From known 
analytic properties of M and U, we note that tp\i^) ^.nd ipxix) are entire in A. These fundamental 
solutions satisfy the Wronskian relation W^'^(a;) = 'W\s{x) with 

m ^ y + (24) 



The corresponding Green's function in equation ( 12 ) is meromorphic with simple poles given by 



the simple zeros of Wy. A — — A„ — — |7i|n, n — 1,2,3, The zeros {A„}„>i can be used as 

initial guesses for finding the eigenvalues {A„}„>i. In particular, for (or p^'~) these zeros are 
initial guesses for the zeros {A„}n>i (w.r.t. A) of W^'^. (oi' ^-\-x+ai^))- To compute 



the Wronskians within the spectral expansion formula (17) we require formulae for the derivatives 



of the fundamental functions in (231 with respect to x. Such derivatives are obtained by using the 
differential recurrences: ^A^(a, b, z) — {a/b)M.{a + l, 6+1, z) and j^U{a, b, z) — —aU (a+1, 6+ 1, z). 
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3.2 Alternative Diffusion Models 

The diffusion canonical transformation is defined as a combination of a change of measure and a non- 
hnear mapping. Consider a solvable underlying time-homogeneous diffusion^ say X^''^ = {xf*^}t>Oi 
defined on the state space T — {l,r), — oo < I < r < oo, and specified by smooth drift, ao{x), 
and diffusion, b{x), coefficients and two fundamental solutions Tpf^ = and (/j^"'' = $a of 

:~ + ao{x)f'{x) — \f{x), A G C, a; G X, subject to appropriate boundary 

conditions. and $p and are respectively increasing and decreasing positive functions of a; S X 
for real values of the parameter p > 0. By applying the change of measure, the solvable underlying 
diffusion is transformed into another diffusion proces^X — {Xt G I}t>o with generator as in (10): 



igf){x) ^ lb'ix)f"ix) + (^aoix) + b'ix)"^^ fix) , (25) 

^ V ' 

—a{x) 

where Up{x) = qi'^p{x) + q2^p{x) is a strictly positive function with constants qi^2 > 0, gi + 52 > 0. 
For a more detailed general discussion of this transformation to various solvable diffusions we refer to 
[U [S] . We simply note here that this corresponds to a (time- homogeneous) Doob-/i transform which 
is generated by the so-called (p-excessive) generating function which is here given by h(x) = Up{x). 

A transition density p of the diffusion X relates to a transition density po of an underlying 
diffusion X'") as follows: 

p{t-x,y) = e-P^'^^Poit-x.v), x,yel,t>0. (26) 

Up[x) 

The speed, Tn(a;), and scale, s{x), densities for X are obtained from the the speed, mo(a;), and scale, 
5o{x), densities of the underlying diffusion via m(a;) = mo(x)Wp(a;) and s(a;) = 5o{x)/u'^{x). The 
above relationship follows similarly when additional killing is introduced. In particular, the transition 
PDFs p^ = Pa^, for the diffusions X^'='=, are related to transition PDFs pQ = Pq^ for underlying 

diffusions X(o)^'±, defined by the generators := (G^"^ - at^>i) and gi°^^ := {0^°^ - at^<e), 

as 

p^{t;x,y)=e-P'^p^it;x,y). (27) 

The main important point of this formula is that closed-form transition PDFs, p^ , are obtained 
directly from the closed-form transition PDFs, pQ , for the solvable underlying process. The latter are 

appli ed t o the simpler analytically solvable 
lake the ob^ 



given by the spectral expansion formulae in Proposition [T 
underlying process i.e. in all terms in equations (16L il7\ and (18 1, we make the obvious 



replacements: ^ vl;, $, ^ mo,s ^ So,P^r,A+. ^ W/,*A+a> W^lj+a ^ W^a A+a' ^a+o,a 

3pon 



<;!,A.<A+a ^ <a+..Wa ^ Wt:t{x)/B{x) ^ ^ W^'^ (x) / Bo{x) . Correspondingly the 



eigenvalues {A„}„>i are defined as the set of increasing simple zeros (w.r.t. A) of W_'^^^ -a(^) ^"-"^ 

-A,-A+qI 



pQ and of Po ■ For instance, in case x < £,y < £ we have 



VK*-* - (£) 



^Such a process X and its generator were denoted more explicitly by {X^^'^'loo and Q^P^ in our previous papers 
[4l[6l|5]- To avoid excessive use of notation, and to maintain consistency in this paper, we drop the superscript p and 
write Xj''' = Xt and Q^P^ = Q. Likewise we write rrip = m and Sp = s. The underlying process is denoted here by 
xt'') and has respective speed and scale functions mo and so- 
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where C^'+ :— j^^x+a xi^)\x=-A '^^^^^ eigenvalues {A„}„>i as the set of increasing simple zeros 
(w.r.t. A) ofT4^*,*„._,(^). 



We remark that the formula in ( 27 ) is also readily derived by making use of the fact that a pair 



of fundamental solutions, 



(p) 



ip\ and 



a(p) - 



for process X with generator in ( 25 1 , for any 



given p > 0, is given by ratios of the known fundamental solutions for the underlying diffusion: 

^x+p{x) , . "^x+pix) 

(t'xix) = ^x{x) - 



Up{x) 



Up{x) 



(29) 



Equation (27) then follows from (17) upon applying the basic Wronskian property: W 

W[f,g]l,x) 



I a 

U ' U 



ix) 



^2(^) • The eigenvalues are found by obtaining the zeros {A„}„>i and then adding p to them: 

An = A„ + p, n > 1, i.e. the eigenspectrum is simply shifted by the positive amount p via the Dooh-/i 
transform as seen in (27). 

Finally, the second main step is to obtain a solvable diffusion S = {St = F{Xt),t > 0}, which 
is used here as an asset price model. This process is defined by a strictly monotonic real-valued 
mapping F with F', F" continuous on I. The mapping F admits the general quotient form: 



Up[X) 



(30) 



The infinitesimal generator of the process S is given by 

1 



{G, f){S)^-a\S)nS)+rSf{S), 



(31) 



where S € Xg = (min{F(/+), F(r— )}, max{F(Z+), F(r— )}). The diffusion coefficient (volatility) 
function is 



a{S) = h{x)r{x)\ ^ KxWK^c^^ ^ +c,^p^,]{x)\ ^^^^^^ 



ul{x) 



(32) 



and r is a real constant such that p + r > 0. The parameter r is equal to the risk-free positive interest 
rate. As above, X = F^^ denotes the inverse map. 



3.2.1 The Bessel-K Model with Killing at an Upper Boundary 

Here we specifically consider a 4-parameter Bessel K-family that arises from an underlying (70- 
dimensional) squared Bessel process (SQB), where we shall assume positive parameters p = ^ — 
1 > and V > Q. By applying the Doob transform with Up{x) = x^^^^K^^ {2^2pxlv) to the 
SQB process X^o' e R+ we obtain a diffusion X € 11+ with generator [Q f)[x) = \v'^xf"{x) + 

7o + t^^x ]]"!^^ ") f'{x). The SQB process has speed and scale densities mo(a;) — ^x^ and So(x) = 



For the Bessel-K model the mapping in ( 30 ) is the strictly increasing function 



{2^2{p + r)x/y) 
F(a;) = c— ^ — , ^ ^ , (33) 



where c and p are independently adjustable positive parameters. F (and its inverse X) maps x G 
(0, 00) and S G (0, 00) into one another. The functions 1^ and denote the modified Bessel 
functions (of order p) of the first and second kind, respectively (see P for definitions and properties). 

For such a Bessel K-family of processes on R+, the origin is regular for p e (0, 1) and exit for 
p>\ (i.e. NONOSC) and the point at infinity is 0-NO natural. Hence, to guarantee that we have 
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a process with NONOSC endpoints (i.e. so that we are in Spectral Category I and Proposition [T] 
applies) we introduce additional killing at some upper level h > and consider X = {Xt}t>o G (0, h). 
Hence, the transformation (33) leads to a family of processes {St — f{Xt)}t>o G (0, H), H = F(h). 
The boundary 5* = is exit if /i > 1 or is a regular (specified as killing) boundary if < /i < 1; 
the boundary 5 = H = F(h) is a killing boundary. We note that one way to deal with the process 
on St G 11+ is to consider the limiting case where the upper boundary is progressively increased to 
infinity. As H ^ oo (i.e. h — > oo), the spectral expansions of the transition PDFs, and hence the 
prices of proportional step options, converge to those for the Bessel K-family on R-|_. 

Differentiating (33 1, and applying differential relations zl'^{z) = fiI^{z) + zI^-^-i{z) and zK'^{z) = 
/xif^(z) — zif^+i (zjTgives the volatility function for the Bessel K-family via equation (32 1: 



/^/^ + i(2^2(p+r)x/iy) 



(34) 



where x = X(S') — F ^{S). The fundamental solutions in (29) follow from those for the SQB process 
with the above mentioned boundary conditions at the left and right endpoints / — and r = h; 



^xix)^x-^^'l^{2V2Xx/iy), 

^xix) = a;-^/2[/^(2V2Ah/t^)i^^(2%/2A^/j/) - K^{2V2)^/:y)l^ 



(35) 



These solutions satisfy the Wronskian relation W^'^{x) 

^x 'x+a ' ^x '\+a ' ^x 'x+a required in the spectral expansion formulas are also easily computed by 
applying the above differential relations. To avoid computations of Bessel functions of complex 
argument, we use the following well-known identities: 



i/p(2\/2Ah/z^)so(x). Other Wronskians 



I^,{ix) = i^J^{x), K^{ix)I^,{iy) - I^{ix)K^{iy) = - {J^{x)Y^,{y) - Y^,{x)J^{y)) , 



where and are the (ordinary) Bessel functions (of order /i) of the first and second kind, 
respectively. Firstly, the zeros {A„}„>i of W'^'^_^ (w.r.t. A) can be computed numerically. These 
zeros are all positive and grow quadratically. In particular, A„ = (i^^/8h)j2^^ where jn,fj. are the 
positive simple zeros of the ordinary Bessel function, i.e. J^(jn,p) = 0. The set {A„}„>i is then 
used as initial guess for computing the eigenvalues {A„}„>i which are the simple zeros (w.r.t. A) for 
either Wronskian W^*3^*„ -a(^) '^^ ^*A*-A+a(^)- Combining all quantities gives us and hence 



by (27) 



3.2.2 The Confluent-U Model 

The confluent-U family of diffusions arises by considering the CIR diffusion as the underlying process 
X*^"^ £ R+ with generator Q'^^^ given by (22). Although this family can be defined for a larger set 

270 



of parameters, here we shall assume positive parameters 70, 7i and define v 



1 > 



0, K 



— 271 



The speed and scale densities for the CIR process are mo(x) ~ {K/'^i)x'^e '^^ and 



So{x) — X ^ ^e'^^ . Applying the Doob transform with generating function Up{x) = U{v, + 1, kx) 
to the CIR process gives us a diffusion process X = {Xt{t>o G ^+ with generator {Q f){x) = 

\v^xf"{x) -f ^7o — 71a; -t- ^'^x^^^fj^ f'{x) 1 where p > 0. For the confluent-U family of models, the 
map F is given by the strictly increasing map 



F(x) 



M{v+ ^^,^l+\,nx) 
U {v, + 1, nx) 



(36) 



where c > 0. Differentiating (36 ), while using j^M{a, b, z) — {a/b)Ai[a+l, b+l,z) and j^U{a, b, z) 
—aU{a -|- 1,6 -I- 1,2:), the volatility function in ([32| for the confluent-U family of processes {St ■ 
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f{Xt)} e R+ takes the form 



lj{S) — CKV\fx 



vM{^,^i+ l,Kx)U{v + l,/i + 2, ra) 



(/i + l)ll{v, fl + 1, K,x) 

where x — X(S') = F^^(S'). For the confluent-U model the origin S* = is exit if /i > 1 and regular 



killing if < /i < 1; the point at infinity is natural. The fundamental solutions in (29 1 follow from 
those for the CIR process: 

■^x{x) = M{\hl,^l+l,Kx), ^x{x) =lliX/ji,fi + l,Kx). (38) 

For the confluent-U model, both endpoints of the state space (0, oo) are NONOSC so that Spectral 
Category I holds and Proposition [l] applies. The functions in (38) satisfy the Wronskian relation 

W^'^';_^ix)/so{x) = with 

r(-A/7i) 

The zeros {A„}„>i (w.r.t. A) are simply given by the poles of the Gamma function in the de- 
nominator. Hence, the eigenvalues grow linearly and are given by A„ = 7i(n — l),n = 1,2,.... 
The set {A„}„>i is then used as initial guess for computing the eigenvalues {A„}„>i which are 
the simple zeros (w.r.t. A) for either Wronskian W^-^_^^ _^{e) or W^-^_^^^{£). The Wronskians 

^x 'x+a' ^x 'x+a' ^x x+a required in the spectral expansion formulas are computed by applying the 
above differential relatio ns f or the confluent M and U functions. Combining all quantities gives us 



= w \ ■ (39) 



Po.a and finally by (27 1. 



3.2.3 The UOU Model 

We now consider the regular Ornstein-Uhlenbeck (OU) process — oo, oo) with constant 

diffusion coefficient b{x) — v > and linear drift coefficient ao{x) = —^ix, 71 > 0. The fundamental 
solutions for this OU process are 

<fxix) = e^-"/^D_x/-,A^V^) and ^^^ix) = $a(-x) = e'^^" /^D_x/^^{-xVk) 

where Di,{-) is Whittaker's parabolic cylinder function (see [T]). 

We now apply the diffusion canonical transformation to the OU process with choice of parameters 
qi = 1,(72 = 0, i.e. with generating function Up(x) = ^p(a;) = e'^^ ^"^D^p/j-^^x^/k) to obtain the 

process X - {X* e R}t>o having the generator (G f)ix) = ^i^^fix) + (-lix + iy^l^) f{x). 



The choice of function in ( 30 1 given by 



F(^) = c-l^f^ = c— ^- , 00, (40) 
maps a; G R onto S G (0, 00) and is monotonically increasing. This transformation gives us a family of 



asset price processes with generator in (31 1 that is referred to as the unbounded Ornstein-Uhlenbeck 



(UOU) model with the diffusion coefficient function given by 

a{S) = c^y^ ' + „ " \ , (41) 

71 D_JL{X^/K) 71 D_J^{Xy/K) D_JL{Xy/K) 

K 71 71 71 J 
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where x = X(5) = F"i(S'). Both ondpoints, S" = and 5 = oo, of the UOU process are NONOSC 
natural boundaries. Hence, Spectral Category I holds and PropositionJ^applies. 

The fundamental functions above satisfy the Wronskian relation W _\ _-J^x) j Zf^{x) — yvi"| where 

(0) ^ ^ 

' r(-A/7i) ^ ^ 

has simple zeros {A„}„>i given by A„ = 7i("- ~ 1)7 = Ij 2, . . .. The set {A„}„>i is then used 
as initial guess for computing the eigenvalues {A„}„>i which are the simple zeros (w.r.t. A) for 
either Wronskian W^x+o^-x^^) or W^-a*_a+q(^)- The other Wronskians W'*^!^^ ' ^*A+a ' ^*A+a 
computed by applying the differential relation ^I3 _^(z) = — (z/2)D_^,(z) — i;L'_,j_i(z). Combining 



all quantities gives us pg'^ and hence p^'* via (27) 



4 Numerical Evaluation of Step Option Prices 

4.1 Spectral Series Expansions 

In this section, we discuss the computational details for computing step option prices using (|3]) and 
the spectral expansion formula for the transition PDFs p = ff^{T\ x, y). Given a discount factor a, 
a level L and spot S (hence £ — X(L) and x = X(5')), the PDF in ^ is to be computed for varying 



values of y G [j/min, J/max]; P IS Computed on such values by truncating the spectral expansion in ( 16 ) 
to the first N terms, where N is sufficiently large: 

N 

pt^iT; X, y) « m(y) ^ e-^"^0l'± (x)0f-± (y). (43) 

n=l 

The numerical procedure consists of the following basic steps. First, compute numerically the 
first N eigenvalues, {A„}i<„<jv associated to the process X. For the solvable asset price processes 
arising from the Doob transform we simply compute the first N eigenvalues {A„}i<„<jv associated 
to the underlying process X^"^ and thereby obtain A„ ~ A„ + p. Note that the eigenvalues grow 
linearly (for the CEV, confluent-U, and UOU models) or quadratically (for the Bessel-K model with 
killing) as n increases. The computations of the terms in the spectral series expansions can be 
split in three parts. That is, we can individually compute the parts that only depend on y, x, and 
£, respectively. Partial derivatives of Wronskians with respect to A can be calculated numerically 
by using a central finite difference approximation. Upon completing the evaluation of the spectral 
series, a quadrature rule (e.g. the adaptive Simpson rule) is applied to compute the integral in (|3|. 
Clearly, this computational scheme can be easily parallelized at different stages. 

For most of the above models, including the CEV, the computation of the spectral expansion for 
the transition PDF p requires many evaluations of the confluent hypergeometric function U{a,b,z) 
for negative values of a. To avoid introducing numerical errors while computing U, the rescaled 
version of the confluent hypergeometric function is used: 

U{a + n,b,z) sin(7r(6-a)) (-1)" T{b-a-n) 



r(-a) sin(7r5) B{~a,b) T{b - a) 



-M{a + n, b, z) 



sin(7ra) ,r(l — n — a) t , , , 

+ -1 " — 6 - 1)^— -^z^-''M{a - 6 + 1 + n, 2 - 6, z , 

TT 1 (—a) 

where a < 0, a ^ Z, and n = 0, 1, 2, . . .. Here B denotes the Beta function. The parabolic cylinder 
function D,^ can be expressed in terms of the Kummer confluent hypergeometric function M provided 
that v ^ X. The rescaled version of D is as follows: 

D.+n{z) ^ „/2 -.V4 r- ( ■M(-(^ + n)/2,l/2,zV2) _ 72zX((l - ^ - n)/2, 3/2, zV2) \ 
r(i./2)-2-/2 [ r{i^/2)r{{l-i^-n)/2) T{iy/2)Ti-{i^ + n)/2) 1' 
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where v > 0, v ^ and n G Z. To compute the products of gamma functions (for large values of 
ly) in the aforementioned formula, we use Euler's reflection formula T(l — z)T(z) = — ^ and the 

^ ' \ / \ / Sin[7TZ) 

following asymptotic series: 

r(x+i) ^,^^11, 5 21 



= x/x 1- 



T{x) V 8a; I28x^ 1024x^ 32768x4 

Since numerical errors can be introduced when directly computing ^^id respective Wron- 

skians for large negative values of A, we can define a rescaled version of (f>\ denoted by 4>x that has 
a better asymptotic behaviour as A — > — oo. For the CEV model we define 

For the Confluent -Z^/ model and UOU model, we respectively set 

<Px(x) — Hi— and (p\(x) = — — A < -p. 

^ 71 / V 271 ' 

4.2 Monte Carlo Bridge Approximation 

We can compare numerical values obtained by using the analytical spectral expansions of the previous 
section with Monte Carlo approximation values. In |16j . a novel algorithm for the exact simulation 
of occupation times for a Brownian bridge is constructed. The method is used to approximately 
sample occupation times for a nonlinear solvable diffusion that admits an exact path simulation. 
Such an approximation sampling algorithm can be applied to the CEV model and other solvable 
hypergeometric diffusions (i.e. the Bessel-K, confluent-U, and UOU models) considered in this paper 
for which an exact path simulation algorithm is available (see [13 )• For example, consider the CEV 
asset price process S. There exists a strictly increasing mapping X that maps S into the CIR diffusion 
X whose volatility is a square-root function, u^/x. The increasing mapping y{x) — -^/x reduces 
the process X to a diffusion Yt ~ Y(Xt) whose diffusion coefficient is equal to one. Thus, we have 

Yt = Y(X(S't)) or St = F(Y-i(i^t)) = ^i'r^?) for ah t > 0, where F = X-\ 

On short time intervals [^1,^2] such a diffusion pinned at points yi and 1/2 at respective times ti 
and ^2 can be approximated by a Brownian bridge from j/i to 1/2 over [ti, t2\. Therefore, occupation 
times of the S bridge process on short time intervals can be well approximated by Brownian bridge 
occupation times. Again, we use the fact that a monotone transformation of a random process does 
not change the occupations times: A^'^ = ^fx^g)^- 

Our approach for the approximate sampling of occupation times A^'^ works as follows. 

1. By using an algorithm from |17j . draw a sample path 6*4 j, . . . , St^,, for a given time partition 

{{) = h<t2<...<tM = T}. 

2. Obtain the respective sample path of the underlying process with unit diffusion coefficient by 
using the transformation Yt- — Y(X(S'(.)) for each i = 0, 1, . . . , M . 

3. Sample the occupation times of A^^^ ^ ^ j for the Brownian bridge from Yt^_-^ to Yt- over ti] 
for each i = 1, . . . , M. Here, £ = Y(X(L)). 



Obtain the approximations A^'^ w Si^i A^'^ 



U-i.ti 



Note that the algorithms developed in jTTl allow us to simultaneously sample the first hitting time 
at zero, tq, and a sample path. If tq < T, then the option is worthless. So the simulation of the 
occupation time can be skipped whenever tq < T . 
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4.3 Numerical Results 

Let us compute prices and deltas of proportional step-down call and step-down put options under 
four different asset price models. The call and put payoff functions are respectively e~""^r' (5't — 
K)+tr„>T and e^^^T (^K ~ S't)+1to>t with proportionality factor a = 5, level L — 90, and 
maturity time T — \. 

Let us begin with the CEV model. It is known from [11] that the negative elasticity values /3 
are typical for stock index options such as S&P 500. We use /3 — —2. The value of b is selected 
so that the local (instantaneous) volatility = (j(Sq)ISq at the spot value Sq — 100 equals 25%: 
(5 = oqS'q^ = 0.25 • 100"-^. For /3 = -2, we have 8 = 2500. The CEV model is compared 
with the other hypergeometric diffusion models with state-dependent volatility function cr(iS'). The 
parameters of the models are adjusted so that the local volatility cr(S'o) at spot 5*0 = 100 is fixed 
at 25%. The parameters of all four asset price models are summarized in Table [T] It should 
be clear that the set used in Table [T] is not the only choice giving (t(S'o) = 25%. In fact, there 
is a continuum of parameter sets for which we can have a fixed value for the local volatility. The 
different parameter values allow us to adjust the steepness, skewness or smile features afforded by the 
various models. This is one important attractive feature of these models, particularly for calibration 



purposes. Figure la illustrates the variety of typical shapes of the local volatility functions (j{S)/S 
when choosing one model over another for given choice of parameters in Table [T] Within a given 
model, we can also further adjust the shapes by varying the model parameters. 
Test 1. 

First, wc calculated the values of step-down call and put options with fixed spot = 100 and varied 
strike K e {80,90,100,110,120} under the four models (as given in Table [T]). The computations 
were done in Matlab with the use the QUADV routine, which numerically evaluates integrals using 
the recursive adaptive Simpson quadrature rule. The absolute error tolerance was set at 10~®. 



Computation of the integrals in (19) and (20 1 with several strikes K reduces to the numerical 



evaluation of integrals of the following two forms 



f F(y) pi'± (T; x,y)dy and / (T; x, y) dy 



on several intervals (a;^, Xi^i). The results of the numerical tests are presented in Table [2] Figure lb 
provides typical shapes of the stock price PDF p{S) :— p^^^ {T; x ,X{S))X' (S) (where T = ^, x = 
X{Sq = 100), and (. = X(L = 90)) computed for the four asset price models. 

All computations were done on a Hewlett-Packard(R) Notebook PC with a four-core Intel(R) 
Core(TM) 17 CPU Q720 @ 1.6GHz and 4 GB of memory. Some details regarding the computational 
time are provided in Table [3j Here we use the following notation: N is the number of terms of the 
truncated spectral expansion, rgp.c xp. is the time required to compute the spectrum {A„}i<„<Ar and 



evaluate the functions 4>n~a{^) in (16), and Tquad. is the time required to numerically evaluate the 
integrals in (19)-(20). Note that me computational time for the Bessel-K model is much smaller 
than in the other models thanks to fast and robust numerical Matlab routines for computing Bessel 
functions. Computations could be drastically sped up if the Matlab routines were translated into 
the machine code, faster and more robust routines for computing confluent hypergeometric functions 
were available, and the code was further optimized. 
Test 2. 

We note that the spectral expansion algorithm allows us to efficiently and simultaneously compute 
option values and deltas for several different strikes and spots. Thus for the second numerical test we 
calculated option values and deltas for a range of spot values. To speed up the computations, we used 
the Simpson quadrature rule with a uniform grid. It allowed us to parallelize the computations of 
spectral expansions by computing individually the parts that only depend on y, x, and respectively. 
The model parameters and problem parameters used here were the same as those for the first test. 
Figures [2]-[5] demonstrate the results obtained. 
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Test 3. 

The results obtained for the CEV model are compared with the Monte Carlo (biased) estimates with 
M = 10^ sample paths and At = 0.05. The results of the numerical tests are presented in Table |4] 
We observe good agreement between the results provided by the spectral expansion method and the 
Monte Carlo algorithm. 
Test 4. 

In another numerical test, we study the sensitivity of the step option price under the Bessel-K 
model as the local volatility function changes its steepness. The steepness of the local volatility was 
controlled by varying the parameter from 0.1 to 0.9. The parameter c was calibrated so that the 
local volatility at 5*0 — 100 is fixed to 25% in all cases. Figure [6] contains both plots of local volatility 
functions and the step-down put option prices. We observe that an increase in the steepness of the 
local volatility tends to decrease the step-down put values. This seems rather intuitive as an increase 
in steepness tends to increase the occupation time below a given level L. 
Test 5. 

Step option prices converge to a barrier option price as a — > oo. In fact, the Green's functions, and 
hence the respective transition PDFs converge to the respective functions for the process having the 
given upper or lower killing barrier level L. In the next numerical example (Figure [7] and Table [5]) 
we show, by computational implementation of the spectral expansions for the transition PDF, how 
the price and delta sensitivity of a step-down call option changes as a increases. The computations 
were done for the Bessel-K model using the parameters in Table [T] The spot price and strike are 
fixed at 100. 
Test 6. 

We studied the convergence of the spectral expansion method as N , the number of terms, increases. 



Figure 8a illustrates the convergence of the PDF given by ( 43 ) as iV increases. The computations 
were done for the Bessel-K model whose parameters are specified in Table [l] The accompanying 
table in Figure [8b| contains the step-down call and put option prices for S[) = K = 100, corresponding 
to using the truncated spectral expansion in (43) for relatively small N number of terms. We hence 
observe typical rapid convergence in the computed prices with the use of the truncated spectral 
expansion formula. 
Test 7. 

Recall that the eigenvalues {A„}„>i grow linearly (for the CEV, confluent-U, and UOU models) or 
quadratically (for the Bessel-K model with killing) as n increases. Therefore, the introduction of a 
killing upper barrier in a hypergeometric diffusion model allows us to accelerate the convergence of 
spectral series expansions. Since the step options are here defined such that they become worthless if 
the upper level h is hit before the maturity time, the option values are biased. In this last numerical 
test we study how such a bias depends on the level h. Figure [9a| contains the graph of the initial 
price Cg^gp(S'o = 100, T = 0.5, — 100) of a step-down call option plotted as a function of the 
level h. As h oo, the probability of hitting the level h decreases and the option price increases, 
asymptotically approaching the option price for the model without killing at an upper level. Again, 
the computations were done for the Bessel-K model. Table 9b gives the put and call option values 
when the upper killing level h changes from 150 to 400. 



5 Conclusions 

One main contribution of this paper is the development of new analytically closed-form spectral 
expansion formulae for the transition probability density function under the CEV model, and under 
various other solvable families of multi-parameter diffusion processes having nonlinear local volatility, 
in the presence of killing at an exponential stopping time (independent of the process) of occupation 
above or below any fixed level. The spectral expansions in Proposition [l] are applicable to a general 
class of diffusions. This paper has successfully implemented the spectral expansions under the CEV 
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model and three other main families of nonlinear local volatility models. As shown in recent papers, 
the nonlinear local volatility models are useful for describing asset price dynamics and for pricing 
standard, lookback and barrier options in finance. This paper further succeeds in providing an 
analytical framework for the risk-neutral pricing of classes of occupation-time options under these 
models. In particular, numerical test results show that the spectral expansions are rapidly convergent 
and provide an efficient method to compute the prices of any proportional step-up and step-down 
options. Moreover, the option Greeks (e.g. delta sensitivity) are also readily and simultaneously 
calculated by simply taking analytical derivatives of the spectral expansions without any loss of 
precision. The spectral expansions converge more rapidly with increasing time and hence the prices 
and Greeks are computed even more efficiently for longer dated options. This offers a significant 
computational advantage in comparison to any Monte Carlo method. 

This paper also derives a general analytical expression for the resolvent kernel (i.e. Green's func- 
tion) of solvable diffusions with killing at an exponential stopping time. This result, by itself, is also 
useful for analytically computing the Laplace transform of certain conditional expectations involving 
functionals of the occupation time for various families of diffusion process above or below a given 
level. In particular, Lemma[T]gives analytical formulae for the conditional expectations in equations 
Q and ([5]) for any solvable diffusion model. For example, these formulae automatically generate 
the expressions tabulated in [2| for various drifted Brownian motions, geometric Brownian motion, 
the (squared) Bessel process, the Ornstein-Uhlenbeck (OU) and radial OU processes. Moreover, by 
the Doob transform employed in this paper, the respective expressions for the Laplace transform 
of the conditional expectations now also extend readily to various newly solvable diffusions. The 
results follow simply from the known fundamental solutions for the underlying process and their 
Wronskians. 

The theoretical development in this paper also sets the foundation for further analytical ex- 
tensions and applications involving the occupation time of newly solvable diffusion processes. For 
example, Lemma [T] can also be extended to cover expectation formula for the case of killing in 
proportion to the occupation time between two levels or a linear combination of occupation times 
below and above a given level. In turn. Proposition [l] can then be extended to cover such cases. 
As long as we are in Spectral Category I, the spectral expansion formulae for the relevant transition 
probability density functions will have a series representation. The inclusion of additionally imposed 
killing, i.e. the usual restrictions on the supremum or infimum of the process by specifying one or 
two interior levels, is also readily handled via the fundamental solutions with appropriately posed 
boundary conditions. 



A Proofs 



A.l The proof of Lemma [T] 

The expectations in Q and ([5| are respectively given by 

d 



dy 



E4e-"^-xi^^^^] ay = AGi'±(x,y,A)dy. 



Hence, by a standard application of the Feynman-Kac formula [e.g. see pages 105-106 in jjj, but here 
generalized to the diffusion with generator defined in ( 10 ) above] the respective Green's functions 
satisfy the ordinary differential equations 

(^ - (A + al,<,))Gi'-(x, y, A) = -Six - y) and (^ - (A + at,>,))&^+ [x, y, A) = -S{x - y), 

where 5{-) is the Dirac delta function. These two equations are solved in the same manner as follows. 

Consider the equation in G^'^. For x < the solution is a linear combination of the pair of 
fundamental solutions {(/'A+a, V-'a+q} of {Q — {\ + oi))Lp = 0. For x > the solution is a linear 
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combination of the pair of fundamental solutions {(j)x,'ijj\} of (Q — X)Lp = 0. A pair {ipXjCpx} of 
fundamental solutions for the equation {Q — {X + a\x<ii))G^or — Oj with the same respective left 
and right boundary conditions as the pair {V'a, </'a}, is hence given by 



x<£, 



Aipxix) + B(f>xix), x>i, 



C^Jx+aix) + D(l)x+a[x), X<£, 



(f>\ix), 



x>£. 



(44) 



The constants A, B, C, D are uniquely determined by requiring that these functions are in C^{I), i.e. 
at X = ^: i^xit-) = V^a(^+), ^'x{£-) = and 4>x{i-) = 0a(^+), i'x{t-) - <^'a(^+). The first 



wt 



it) 



set of conditions is a hnear system of equations in A, B with solution A — " ^'Iti^'u^ ■• ^ — M/^°'^rc i 



WT 



and the second set of conditions is a linear system in C, D with solution C = — a+q^ 
w 



Computing the Wronskian for the pair in ( 44 1 gives 



W[^x,i^x]ix) ^ 



S{x) = yvx5{x). 



Here we used the identity W^^':f{x) = W-yS{x) so that VF^>(x)/VF.f;^^(^) = s(x)/s(£) for any 7 e C. 
The Green's function G^'" is then simply given by combining this Wronskian with (44): 



Ga ix,y,X) _ ipx{x hy)(t>x{xy y) 



m(y) 



Wa 



■>Jjx{x A y)(t)xix V y) 



A,A+Q 



(^) 



(45) 



for x,y €X. The Green's function in (45 1 can be recast as in equation ( [15| ) by using (44 1, and the 
definition for G in (12 ), for the four different cases: x,y < £, x < £ < y, x > £ > y or x,y > £. 

The derivation for G„'+ in equation (14) follows the same steps as above. A general solution to 
the equation (t/ — (A + atx>£))Ga'^ = is a linear combination of the pair {V'a, '/'a}, for x < £, and 
of the pair {(j)x+a, V'a+q} for x > £. In analogy with equation (44), a pair of solutions {ipxi 0a}i with 
the same left and right boundary conditions as the pair {iIjx,4>\} is as follows: 



i^xix) 



i^xix), 



X < 



A-ipX+aix) + B(j)x+aix), X > £, 



4'\{x) 



C'ipxix) + D(j>xix), x<£, 



0A+a(a;), 



> £. 



(46) 



The requirement that these functions are in C^(I) then uniquely gives A — 



'^'^ A,A + 



A + a,A + 

computed to be 



, B 



The Wronskian of the two solutions in (46) is readily 



,{£) 



W[^xMx) = 



a{x) = Wx5{x). 



Combining this Wronskian with (46) gives the Green's function: 

G^'+(x,y,A) iJx{x ^y)^x{xy y) ^^^pxix A y)(l)x{x V y) 



m(y) 



Wa 



(47) 
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for y G I. The Green's function takes the more expHcit form in (14) by using (46), and the 



definition for G in ( 12 ), for the four different cases: x,y<i, x<i<y, x>£>yoTx,y>£. This 
completes the proof. 

[Remark: In the limit a oo, the Green's functions G^'^ can be proven to converge to the 
respective Green's functions for the process with killing at an upper (or lower) level £. We do not 
give a proof here, as it is based on the a oo formal asymptotic analysis of the fundamental 
solutions ipx+aix) and 0a+q(^c). By the leading term asymptotics, ^p\-^.a{x)/'lp'^^^{x) — and 

these limits and the asymptotic properties we can arrive at the asymptotic forms for the Green's 
functions in equations (14) and (15) of Lemma [l] In particular, as a — ?■ oo: 



(t>x+aix)/<j)'x+ai^) -> 0. Hence, for example, ^^^^ ^^^^ 



G'^+{x,y,\) ^ G''+{x,y,X) _ 4'x{x Ay)S{xy yJ-X) 
m{y) m{y) WaV-aW 

with G^'+ {x,y, \) = if X > £ OT y > e, and 

Gi--(x,y,A) ^ G''-{x,y,X) _ a; A y; A)0a(x V y) 



m{y) 



x,y <£, 



x,y > £, 



with G^~{x,y,X) = if x < £ or y < £. Here, G^'+(a:, y, A), or G^'~(x,y, A), are the respective 
Green's function for the process X < £, or X > with killing imposed at the upper, or lower, level 
£. The generalized cylinder function is defined by 

S{x,y;X) := (t)x{x)iljx{y) - ■ipx{x)(t)x{y)- 
For real A > 0, S{x, £; A) {S{£, x; A)) is a decreasing (increasing) positive function for x < £ {x > £).] 



A. 2 The proof of Proposition [T] 

Due to the similar structure of the Green's functions Gl'+ and Gi 



as observed in equations ( 14 ) 



and (15), we will only present the proof for p^'^, i.e. equations (16) and (17). The PDF is given by 
the Laplace inverse Pa~^{t; x, y) = ^G^'"''(x, y, A)^ (i) which can be computed for all four separate 
cases in equation ( 14 ) with standard use of the Residue Theorem upon closing the Bromwich contour 



integral on the left-half of the complex A plane. In particular. 



p^'+(t;x,y) 
m(y) 



OO 

^e-^"*Res 



G'^+{x,y,X) 
m(y) 



A - -A,: 



where Res 



m(y) ' - -^'i 



Consider the case where x < £ < y, i.e 



Gi; + ix,y.\} _ 



-ilJx{x)<l>x+a{y)- The only singular- 



ities are the real simple zeros at A 
increasing sequence of eigenvalues solving W"'- 



m{y) M'A+t,A(^)' 

^n^n > 1- of the Wronskian W^^^ a(^)' 'where {A„}„>i is an 



.-A. 



{£) = 0. By analyticity of the fundamental 



functions (w.r.t. A), the residue of the Green's function at these simple poles is then given by: 



Res 



G'^+{x,y,X) 



■ s(£) Res 



■0A(a;)(/'A+a(y) 



:S(£) 



V'_A„(2:)(/>_^ ,„(y) 



W^a^1aWIa=-a„ 
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The case where y < i < x is similar, where 



Res 



A — — A,i 



: s{e) Res 



0A+a(a;)V'A(y) 



A — — A„ 



We now consider the case where x < £,y < £, i.e. 



G'^+{x,y,X) ^ G{x,y,X) ^ Wt:t+M M^)My) 
miy) m{y) W^^^^je) Wa ' 



with 

My) 



given by equation ( jl2| ). Hence, this Green's function has two sets of singularities. 
One is the set of simple poles corresponding to the simple zeros solving Wa=-a„ = 0, i.e. {A,i}„>i 
denotes the eigenvalue set for the Sturm-Liouville problem with generator Q for diffusion X on I. 
The other is the set of zeros A = — A„, n > 1, of W^]^^ a(^)- "^"^^ establish that the only nonzero 
residues are for the set A = — A„,n > 1. Assume the set {A„}„>i is isolated from the set {A„}„>i. 
Then, computing the residue at every simple pole A = — A„ gives: 



Res 



G'^ + {x,y,X) 



Res 



Gix,y,X) . 



Res 



;A 



For A = — A„, the fundamental functions tpx and (f>x are proportional to each other, i.e. (f)^\^{x) = 
j4„^-a„(2;), for some constant An 7^ 0. Hence, the ratio of Wronskians in the above second residue 



term evaluates to 



wtit . W 



Denoting C„ := j^V^xlx^-x ' ^^"^ second residue 



evaluates to 
Res 



X — — A„ 



= -7^V'-A„(a;)V'-A„(2/) = -(/)„(x)(/)„(y) 



where (/)„(a;) :— ±y ^ip_x^{x) is the n-th eigenfunction of —Q for x £l. The first residue has the 
standard eigenfunction product form: 



Res 



Gix,y,X) 

My) ' 



X 



-A,] 



4>n{x)(j)n{y). 



Adding the two terms gives a zero residue at every A = — A„, i.e. Res 



= 0. 



The only nonzero residues are hence due to the assumed simple poles A = — A„, n > 1 and these 



are 



Res 



G'^+{x,y,X} ^ 



Res 



w^rj+affl Mx)My) 



w 



A,A+Q 



(e) 



^-\J^)'>l'-\Jy) 



-A„ 



which is the form in equation (17) for a; < ^, y < ^- We note that if the n-th eigenvalue A„ happens 
to also coincide with an eigenvalue in the set {A„}„>i, say A„ = A™ for some m > 1, then the 
above formula is interpreted as a limit A — > — A„. We remark that for a coalescence of zeros (w.r.t. 
A) of both Wronskians, Wx and a(^)' ^^e point A = — A„ is still a first order pole since 
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(px = Anipx, and hence the Wronskian W^'^_^_^{£) in the numerator is proportional to ^{£) in 

the denominator, at A = — A„. 

The last case where x > i,y > £ follows in very similar fashion. Again, the residues for the set 
A = — A„,n > 1, are all zero and the only nonzero residues are due to the assumed simple poles 
A = — A„, n > 1, where 



Res 



; A — — A„ 



= Res 



r4>,'4> 
A,A+ 



A — —Xn 



-A„+a 



A=-A„ 



Again, the same above remarks apply here as for the previous expression just above. 
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Table 1: Parameters of the four asset price models 



Model Parameters 



CEV 


6 


= 2500; /3 = -2; r = 0.02 






Bessel-K (BK) 




= 0.5; 70 = 2.2; p = 0.00001; c = 


728.7467627; 


h = 500; r = 0.02 


Confluent-U (CU) 




= 0.5; c = 133.1173736; p = 0.01: 


; 1/ = \/2; 71 = 


= 0.1; r = 0.02 


UOU 


P 


= 0.001; z/ = 2; c = 71.11606167; 


71 = 0.2; r = 


0.02 



Table 2: Values of the step-down call and put options computed for a range of strikes under the four 
asset price models. The parameters used are = 100, T = 0.5, a = 5, i = 90. 



Step Calls 

K CEV BK CU UOU 

80 20.364424 19.774476 20.657331 20.603226 

90 13.359199 12.953074 13.563118 13.043498 

100 7.336247 7.351327 7.403872 7.244763 

110 3.130114 3.610598 3.107478 3.844225 

120 0.948158 1.548517 0.947478 2.057045 



Step Puts 



CEV 

0.295586 
0.840873 
2.368432 
5.712811 



BK 

0.167632 
0.748891 
2.549805 
6.211737 



CU 

0.340048 
0.859860 
2.314636 
5.632266 



UOU 

0.013590 
0.379618 
2.506640 
7.031858 



11.081367 11.552317 11.086289 13.170435 



Table 3: Computational costs of numerical evaluation of step-down call and put options (for one 
spot value and five strike prices) under the CEV, CU, UOU, and BK asset price modc^ls. The Matlab 
code was run on a Hewlett-Packard(R) Notebook PC with a four-core Intel(R) Corc(TM) i7 CPU 
Q720 @ 1.6GHz and 4 GB of memory. 



Model 


TV 


T 

sp.exp. 


-^quad. 


CEV 


150 


8.81 sec 


60.70 sec 


BK 


50 


0.77 sec 


2.94 sec 


CU 


300 


20.96 sec 


148.68 sec 


UOU 


150 


5.36 sec 


52.63 sec 



Table 4: The Monte Carlo biased estimates of proportional step-down call and put prices under the 
CEV model using the Brownian bridge interpolation method are tabulated for various values of strike 
K. Monte Carlo estimates are compared with analytical approximations provided by the spectral 
expansion method. Here, sm denotes the stochastic error. The parameters used are 5*0 = 100, 
r = 1, r = 0.1, 5 = 2.5, /3 = -0.5, a = 0.5, L = 90, Ai = 0.05. The number of sample paths is 
M = 10^ 



K 


Step Calls 

MCM Estimate ± Sm Analyt. Approx. 


Step Puts 

MCM Estimate ±sm Analyt. Approx. 


90 
100 
110 


20.9950 ±0.0009 20.993325 

14.8192 ±0.0006 14.817208 
9.8621 ±0.0004 9.860234 


2.0416 ±0.0006 2.039807 

4.1988 ±0.0010 4.195828 
7.5750 ±0.0017 7.570991 
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Table 5: Step-down call and put option values are computed under the Bessel-K model for increasing 
values of a. The option parameters are K — 100, T — 0.5, L = 90. The model parameters are 
specified in Table [T] The case with a = oo corresponds to the double knock-out barrier option with 
barriers L = 90 and U = 400. When a = 0, the step call and put options reduce to the European 
call and put options, respectively. 



a 


Call Value 


Put Value 





7.525593 


6.530576 


1 


7.483054 


5.213809 


5 


7.351327 


2.549805 


10 


7.240869 


1.469595 


25 


7.060945 


0.752496 


50 


6.925459 


0.524287 


100 


6.806920 


0.404356 


200 


6.710443 


0.336245 


500 


6.615400 


0.285627 


1000 


6.563974 


0.263218 


oo 


6.494245 


0.218820 




(a) The local volatility a(S)/S. (b) The PDF pi'". 

Figure 1: The local volatility function a{S)/S and corresponding transition PDFs p^^, as function 
of spot S, for the asset price process with killing at an exponential stopping time of occupation 
below a fixed level L, are computed for four asset price models specified in Table [l] The PDFs p^^~ 
are computed for the following parameters: Sq = 100, a = 5, L = 90, and T = ^. 



G. Campolieti, R. Makarov, and K. Wouterloot, Pricing Step Options 



27 




(a) Put values under the CEV model. (b) Put values under the UOU model. 



Figure 2: Values of the step-down put option, as function of spot S, computed for a range of strikes 
under the CEV and UOU models. 




(a) Call values under the CU model. (b) Call values under the BK model. 



Figure 3: Values of the step-down call option, as function of spot S, computed for a range of strikes 
under the confluent-U and Bessel-K models. 
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(a) Put deltas under the CEV model. 



P 0.2 




(b) Put deltas under the UOU model. 



Figure 4: Deltas of the step-down put option, as function of spot S, computed for a range of strikes 
under the CEV and UOU models. 




(a) Call deltas under the CU model. (b) Call deltas under the BK model. 



Figure 5: Deltas of the step-down call option, as function of spot S, computed for a range of strikes 
under the (a) confluent-U and (b) Bessel-K models. 
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(a) Local volatility functions. 





= 0.1 




= 0.3 




= 0.5 


n 


= 0.7 




= 0.9 






(b) Step put values. 



Figure 6: Local volatility functions and put option prices computed, as function of spot S, under the 
Bessel-K model for a range of values of fj,. The parameter c varies accordingly so that (t(100)/100 = 
0.25 for all choices of fj,. The other model parameters are as specified in Table [T] The other 
parameters are T — ^, a = 5, L ^ 90, and K ~ Sq = 100. 



Values of the Step Call for Bessel-K 




80 85 90 95 100 80 85 90 95 100 

S S 



(a) Call values. (b) Call deltas. 

Figure 7: Step-down call option values and deltas are computed under the Bessel-K model (as 
specified in Table [T]) for increasing values of a. The option parameters are K = 100, T = 0.5, 
L = 90. 
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N = 


2 








4 






N = 


6 


0.02 






10 




150 
S 



N 


Call Value 


Put Value 


4 


0.818540 


1.705722 


6 


6.290608 


2.867394 


10 


7.593736 


2.745855 


15 


7.377383 


2.539072 


20 


7.351545 


2.550244 


50 


7.351327 


2.549805 



(a) Approximations of the PDF 



(b) Call and put values. 



Figure 8: The convergence of truncated series approximations of the PDF p^^", as the number of 
terms N increases. The computations were done for the Bessel-K model whose parameters are 
specified in Table [ij The other parameters are T = i, a = 5, L = 90, and S'o = 100. The step-down 
option values are computed for the strike price K = 100. 



h Call Value Put Value 




6.594098 2.549803 

7.076232 2.549805 

7.261562 2.549805 

7.324466 2.549805 

7.343837 2.549805 

7.349354 2.549805 

7.351269 2.549805 

7.351326 2.549805 

7.351327 2.549805 
7.351327 2.549805 
7.351327 2.549805 
7.351327 2.549805 



150 200 250 300 350 400 



Upper barrier h 

(a) The call value as function of h. (b) Call and put prices. 

Figure 9: Convergence of the prices C^gp{So — 100, T = 0.5, K — 100) of the step-down call 
option as the imposed killing level h increases. The computations are for the Bessel-K model whose 
parameters are specified in Table [T| The other parameters are a — 5 and L = 90. 



